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We calculate the odd-parity, radiative {i > 2) parts of the metric perturbation in Lorenz gauge 
caused by a small compact object in eccentric orbit about a Schwarzschild black hole. The Lorenz 
gauge solution is found via gauge transformation from a corresponding one in Regge- Wheeler gauge. 
Like the Regge- Wheeler gauge solution itself, the gauge generator is computed in the frequency 
domain and transferred to the time domain. The wave equation for the gauge generator has a 
source with a compact, moving delta-function term and a discontinuous non-compact term. The 
former term allows the method of extended homogeneous solutions to be applied (which circumvents 
the Gibbs phenomenon). The latter has required the development of new means to use frequency 
domain methods and yet be able to transfer to the time domain while avoiding Gibbs problems. 
Two new methods are developed to achieve this: a partial annihilator method and a method of 
extended particular solutions. We detail these methods and show their application in calculating 
the odd-parity gauge generator and Lorenz gauge metric perturbations. A subsequent paper will 
apply these methods to the harder task of computing the even-parity parts of the gauge generator. 
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I. INTRODUCTION 

The increasing maturity of space-based gravitational wave detector concepts [HH] has in part motivated considerable 
work in the last fifteen years on self-consistent calculations of extreme-mass-ratio inspirals (EMRIs). Such a system 
consists of a small compact object of mass ^ ~ 1 — IOMq (e.g., neutron star or black hole) moving on a decaying 
orbit about, and ultimately into, a supermassive black hole of mass M ~ 10^ — IO^Mq ^ ^. Irrespective of when a 
detector might launch, there is also simply considerable theoretical interest in the problem of motion of a point mass 
in a background geometry in general relativity, influenced by its own self- force [3]. 

The extreme mass-ratio lends itself to use of black hole perturbation theory. In the limit of /i the small mass 
orbits on a geodesic of the massive black hole background with constants of motion. At next order the small mass 
draws up a small perturbation in the metric, which results in gravitational radiation fluxing to infinity and down the 
horizon of the massive black hole. The metric perturbation (MP) also acts back on the small body locally (through 
a self- force), giving rise to dissipative effects that cause the orbit to decay and to small conservative corrections to 
the motion. The perturbation problem is singular in several respects [4j, with a divergence in the MP at the particle 
location making the motion correction also formally divergent. A general understanding of how to treat the self- force 
(i.e., regularize it) in an arbitrary spacetime was given by Mino, Sasaki, and Tanaka [5] and Quinn and Wald [5]. 
Practical procedures for regularizing the self- force in numerical calculations followed (e.g., [7j). 

The physical retarded field p™* can be conveniently split into regular (i?) and singular (S) parts, as first introduced 
by Detweiler and Whiting [8]. The advantage of this split is that while the singular contribution to the MP pj^j, 
satisfies the inhomogeneous field equations, it does not contribute at all to the self-force. On the other hand, the 
regular contribution, p^^, is a smooth, homogeneous solution to the field equations and, through a projected gradient, 
is entirely responsible for the self-force. Indeed, when interpreted this way, the regular field can be thought of as an 
external field which sources the deviation from geodesic motion on the background metric g^^. The motion of the 
particle is then geodesic on the spacetime The singular part of the MP is calculated analytically in Lorenz 

gauge, and an expansion provides the regularization parameters [7] . The singular part is then subtracted from the full 
retarded field mode by mode in a spherical harmonic expansion, allowing the difference {p^^) to converge. While in 
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principle [51 [TO] the full retarded field could be calculated in a variety of gauges, in practice most calculations [TTHTOj 
have also used Lorenz gauge to find p^^^ . 

We are developing a set of techniques and assembling a computer code to calculate with high accuracy the first-order 
MPs from a small compact object in a generic orbit about a Schwarzschild black hole. The need for high accuracy is 
related to a set of arguments that have been made for years that EMRIs should be calculated through second-order 
in perturbation theory |3l IMUTT] . One particular argument centers on calculating the phase evolution of an EMRI 
and using it in interpreting data from a detector. For a small mass ratio e = ^/M (and in the absence of transient 
resonances [TB] which may occur during EMRI evolution on a Kerr black hole), we expect that as an EMRI evolves 
through a detector passband the gravitational waveform will accumulate (schematically) a phase of 

$ = Ki- + K2e" + ^36^ + • • • , (1.1) 
e 

where the k's are coefficients of order unity that depend upon, among other things, the lower and upper limits in 
frequency of the detector response. The first term refiects the dissipative effects of the first-order self-force in spurring 
a decay of the orbit. The second term would result from second order in perturbation theory. For example, with 
€ ~ 10~^, an EMRI might be observed to accumulate a total phase of $ ~ 10^. For matched filtering purposes we 
might need to compute the phase to an accuracy (5$ < 0.1, and thus a fractional error of < 10^^. However, the error 
in phase in using the first-order calculation alone is ~ Hence, the need for a second approximation to take full 

advantage of a detector output. But there is a corollary to this argument. We cannot possibly hope to make use of a 
second-order calculation if we have not already computed the first-order self-force to a relative accuracy much better 
than 0{e). The requirement might be at least several orders of magnitude better than 0{e) to make a second order 
calculation worthwhile (say 10~^ to 10~* in the example). Furthermore, given that computation of the self- force is a 
numerically subtractive procedure, the first-order pre-regularization field contributions likely need to be known even 
more accurately (perhaps 10""'^^ to 10~^°). 

An accurate calculation strongly suggests use of Fourier decomposition and frequency domain (FD) methods, to 
gain the benefit of integrating ordinary differential equations. Ultimately we are interested in the time-dependent 
self-force and must transfer back to the time domain (TD). For that step, a lynchpin of the effort has been use of the 
recently developed method of extended homogeneous solutions (EHS) [19 , which allows partial Fourier series sums 
for the perturbations to avoid the Gibbs phenomenon and to converge exponentially even at the location of the point 
mass and despite loss of differentiability there. This conclusion has guided others as well |2D1 HI] . 

Like other recent calculations, we want to determine the first-order MP in Lorenz gauge. However, our approach 
is indirect. In an earlier paper |22], we calculated radiative modes {i > 2) of the MP in Regge- Wheeler (RW) 
gauge by applying EHS to solutions of the master equations of the Regge- Wheeler- Zerilli (RWZ) formalism and then 
determining the metric from the master functions. The new aspect of that work was being able to use FD techniques 
and nevertheless determine the metric amplitudes in the TD with accuracy right up to the location of the particle 
r = rp(t), an essential requirement for computing the self- force accurately. EHS works by recognizing that the solution 
of a master equation with a moving singular source (in the RWZ case the source has both a delta function term and a 
derivative of delta function term) is a weak solution of the form ^'(t, r) = ^'"'"(t, r)9[r — rp{t)] + vE'~(i, r) [rp{t) — r], 
where 5'+(i,r) and 4'~(t, r) are differentiable solutions to the source- free master equation. and are in turn 
obtained as Fourier sums of properly-normalized Fourier-harmonic modes that solve the source-free master equation 
in the FD. Since the functions in the separate Fourier sums are smooth everywhere, the lack of differentiability of ^ 
is captured entirely by the 6 functions. We can in turn then calculate the MPs from The result, however, is in RW 
gauge. This paper, and a subsequent one, address the calculation of the infinitesimal gauge generator that transforms 
pJJ^ in RW gauge to its counterpart p^^ in Lorenz gauge. 

This paper is restricted to finding the odd-parity part of the gauge generator, which for each spherical harmonic 
mode has an amplitude that is a solution to a single inhomogeneous wave equation. While the wave equation is simple 
to express, what is more challenging is to find a way to generalize the underlying idea behind EHS to equations with 
non-compact source terms. A substantial part of this paper is devoted to laying out two new analytic/numerical 
methods (method of partial annihilators and method of extended particular solutions (EPS)) we have developed for 
solving differential equations of this type. These techniques will play a central role in a subsequent paper where 
we present the more involved procedure, based on analysis by Sago, Nakano, and Sasaki [23], for determining the 
even-parity parts of the RW-to-Lorenz gauge generator. 

The discussion above begs the question, why the need to transform to Lorenz gauge? In part we know that the MP 
amplitudes in Lorenz gauge are C° in the TD at the particle location. As we discussed in Ref. ^2] the MP amplitudes 
in RW gauge are one or two orders of continuity worse behaved (i.e., some amplitudes are C~^, or discontinuous, and 
some have a delta function term at r = rp{t)). Fig. [l] demonstrates this problem graphically for the £ = 2,m — 1 (odd 
parity) RW amplitudes h^'^ and ft,^'^. The insets show discontinuities in the MP amplitudes at the particle's location. 
The plots also show the cx r growth in the wave pulse amplitude as r — )• 00, refiecting the fact that RW gauge is not 




FIG. 1. The £ = 2, m — 1 mode of the Regge- Wheeler gauge MP amphtudes ht and ■ The orbit is parameterized 
by eccentricity e = 0.764124 and semi-latus rectum p = 8.75455. The plots show the real and imaginary parts of the MP 
amplitudes at time t = 93.58 (where t = is at the periapsis). Dotted vertical lines indicate limits of the source libration 
region. Insets show the discontinuities at the location of the particle. The lack of asymptotic flatness is evident as r — >■ cxd. We 
plot fh^ so that the wave behavior near the horizon can be seen. 



asymptotically flat 24J. Transformation to Lorenz gauge not only improves the behavior of the modes at r = ^p{i), 
it also removes the non-asymptotically-flat behavior seen in the radiative modes of RW gauge. 

Lastly, we note that for non-radiative {I — 0, 1) modes, the RWZ formalism breaks down. Various researchers have 
used different gauges to solve the Einstein equations for these modes. Considering generic motion on a Schwarzschild 
background, Zerilli [2 5) solved for these modes analytically in his own gauge which exhibits behavior for certain 
modes. For circular orbits Detweiler and Poisson [26' showed how to transform Zerilli's solutions to Lorenz gauge. 
Barack and Sago [13] solved the Lorenz gauge field equations directly for these modes (and higher radiative modes). 
The non-radiative modes provide a crucial contribution to the conservative piece of the self-force. In this paper, 
though, we are solely concerned with transforming our RW gauge (£ > 2) MP amplitudes to Lorenz gauge. However, 
we note that while RW gauge is not defined for t <2, there is no such restriction on the gauge transformation. Indeed, 
the work presented here can be directly extended to handle the transformations of non-radiative mode solutions from 
other gauges to Lorenz. 

Throughout this paper we use the sign conventions and notation of Misner, Thorne, and Wheeler |27| and use 
units in which c = G = \. The background geometry is a non-rotating black hole which is described in terms of 
Schwarzschild coordinates — {t,r,d,Lp). 



II. FORMALISM 



Consider the motion of a small compact object of mass /i in orbit about a static black hole of mass Af , with 
/i/M <C 1. The small body gives rise to a perturbation p^^ in the metric relative to the Schwarzschild background. 



ds^ = g^^dx^'dx•' = -f{r)dt^ + f{r)'^dr^ + {d6^ + sin^ Odip^) , 



(2.1) 
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where /(r) 1 — 2M /r. We are concerned here only with the first-order part pl^l of the MP in an expansion in powers 
of /i/M and accordingly simply set p^^, = pjly . The full non-stationary metric is then g^^, — g^^ + p^^{t, r, 0, ip). It 
often proves convenient to work with the trace-reversed MP 

= - Idt^i^Papg"'^- (2.2) 

The Einstein field equations can be expanded about the background geometry. Once the stress-energy tensor 
associated with the motion of the small compact object has been specified, the first-order linear equations 

GS(p-,.) = 87rr(°\ (2.3) 

can be solved to determine p^y The small body is approximated as a point particle and the stress tensor is that of 
the particle moving on a geodesic of the Schwarzschild black hole (i.e., zeroth-ordcr approximation). 

In this work, we are interested in bound eccentric motion. The orbit is parameterized by a pair of constants, which 
can alternately be taken to be Darwin's [5S] eccentricity e and dimensionless semi-latus rectum or the bounding 
values of radial motion Tmin and "rmax, or the specific energy E and angular momentum C. The first integrals of motion 
are integrated using Darwin's curve parameter x s-^d are converted to functions of the form r = rp(t) , ip = fp{t)^ and 
T = Tp{t) (with 9 = tt/2). The orbit has two fundamental frequencies with a rate fir associated with radial libration 
and a mean rate of azimuthal advance fltp. See Darwin [28 for the orbital integration and Cutler et al. [29] for first 
application to black hole perturbations from bound motion. 



A. Harmonic decomposition of the field equations 



The usual route to solving the field equations (2.3) on a Schwarzschild background is to introduce tensor spherical 
harmonics and decompose the equations into individual angular harmonic modes. This approach was followed by 
Regge and Wheeler [3D] and Zerilli [25' in solving for certain master functions that represent the odd- and even-parity 
parts, respectively, of the gravitational field. The MP is then derived from these master functions in Regge- Wheeler 
gauge. In Ref. [22], we used a variant of this approach, with slightly different versions of the two master functions, 
and combined it with a new analytic/numerical method for finding convergent solutions near the location of the point 
particle. 

In that paper we followed Martel and Poisson [31] in using their definitions of the tensor spherical harmonics. We 
recap those definitions here. The unit two-sphere is covered by coordinates {0,ip). Upper-case Latin indices A, 
B, etc. denote these two angular coordinates and associated tensor components. The coordinates {t, r) cover the 
submanifold Ai^. For these coordinates and tensor components we use lower-case Latin indices a, b, etc. The full 
Schwarzschild spacetime is = x S^. The usual scalar spherical harmonic functions are Y^"^{6,Lp). From 
these we define even-parity = D^F^™ and odd-parity = —£a^ DbY''"^ vector harmonics, where Da is 

the covariant derivative on 5^. The metric on the unit sphere is ^ab = diag[l, sin^ 0] and the Levi-Civita tensor 
is Sab- Both are compatible with Da- Dc^ab — DcSab = 0- Martel and Poisson define the two even-parity 
tensor spherical harmonics as QabY^"^ and F^'g = [DaDb + + 1)] y^'", where the latter is trace-free and thus 
differs from that used by Regge and Wheeler {DaDbY^"^) ^30j. They take the odd-parity tensor harmonic to be 
X% = -i {ea^Db + Eb'^Da) DcY'^"', which differs by a minus sign from that of [30]. 

Using these definitions the MP is split into even- and odd-parity parts. Each parity is decomposed into sums over 
£ and m of their respective harmonics. For even-parity there are seven amplitudes that Martel and Poisson define, 
which can be related to those of [3D] and [5S]. They are hu = fHo, htr ~ Hi, hrr = H2/ f, jt — ho, jr = hi, 
Ghoro — Grwj a-nd inhere = ^RW ~ ^(^ + l)G/2. (Note: here and on many occasions later in the text we suppress the 
spherical harmonic scripts £ and m for brevity when no confusion should arise.) Then the even-parity MP is 

Pab^Y.^a'^^'"'' PaB^Y.^i^'^B"'' P^B = ^ (x^^l^^sF^™ + G^^rj-g) . (2.4) 

£,m £,m £,m 

For odd-parity there are three amplitudes, ht = ho, hr = hi, and — —h^^, which are equivalent to the original 
definitions up to sign. The odd-parity MP is then 

PaB = E , PAB = Y ^^^is , (2.5) 

along with the fact that pab — 0. For the balance of this paper we are only concerned with odd-parity. 
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Regge- Wheeler gauge places the (algebraic) condition on the metric that = 0. In this gauge the odd-parity 
field equations become 

2 5^/,RW ^(^ + l)r-4M, 
r 



-dtdrh 



RW 



^2^,,KW^(i+2)(£-l)/,,^ 



r 
1 

'-'t'h 



/ 



dthf^ + fdrk 



RW 



2M, 



(2.6) 



where the source amplitudes P*, P*", and P (for each I and m) are odd-parity projections of the stress tensor, 



P(t,r) = le^r^^^ 



t^^xXb dn. 



(2.7) 



{i + 2y. 

We use an asterisk to denote complex conjugation. The source amplitudes in turn satisfy the contracted Bianchi 
identity 



Qpt Qpr 



ie-m + 2) 



P = 0, 



(2.8) 



dr r 

and the stress tensor itself is taken to be that of a particle in geodesic motion on the background geometry. 

While in principle the coupled equations (2.6) might be solved to yield the metric in RW gauge, the usual approach 
involves defining and using one of several odd-parity master functions and solving a lone wave equation (master 
equation) for this function. The odd-parity part of the metric then is derived from the master function. An equivalent 
master function representation is used for even-parity in RW gauge. In our previous paper [22] we used the odd-parity 
Cunningham-Price-Moncrief (CPM) function, which we refer to here as '^o- This gauge- invariant master function is 
defined in Regge- Wheeler gauge by 



*o(i,r-) 



2r 



^rh^ 



RW 



,RW 



In terms of the tortoise coordinate = r + 2Af ln(r/2M — 1), satisfies the wave equation 

where W2 is the spin-2 Regge- Wheeler operator, a particular case of the spin-s operator 

d"^ . 92 j£(£+l) 2M(l-s2)" 



dt"^ dr'i 



f 



(2.9) 
(2.10) 
(2.11) 



Later, we will have need for the Fourier transform of this operator (dt — >■ — iw), which we will denote Cs- The source 
term for the master equation involves a combination of moments of the stress tensor, 



2rf 



(£-!){£ + 2) 



f 



Go{t) 5 [r - rp(t)] + Fo{t) S' [r - rp{t)] , 



(2.12) 



and is a distribution (see [22] for details) . Once the CPM master function is known the MP amplitudes in RW gauge 
can be reconstructed via the expressions 



f. 



i^r Khr) 2f '"^ {£-!){£ + 2) f ' 



(2.13) 



2-"'-'"' {£-!){£ + 2)' 

Their numerical determination in the time domain with a convergent and accurate behavior everywhere including the 
vicinity of the moving particle was the subject of our previous paper. 

Because we have reason to consider it in what follows, it is worthwhile noting that the original master function of 
Regge and Wheeler, V&rw, is not the CPM master function we use. They are related by 



1 



The RW master function satisfies an almost identical wave equation, 

W2«'Rw(i,'^) = Snw, 

with the only difference being the source term 



pr 



SRwit,r) 



f 



-P^ + fdrP - 



1 - 



3Af 



P 



(2.14) 
(2.15) 
(2.16) 
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B. Gauge transformations 



The exact form of the field equations (2.3) wih depend upon specifying a gauge. As mentioned in the Introduction, 
two frequent choices are Regge- Wheeler (RW) gauge and Lorenz (L) gauge. The small gauge generator S'' that 

+ between the two gauges is on the same order of magnitude as the MP, 



transforms the coordinates 



'-RW 



that is |S^| ~ \p^^\ 1- The MP then transforms as 



Pfiu P^iu ^l'\^l^9^^v^ |aj 



(2.17) 



where stroke |/i indicates covariant differentiation with respect to the background metric. Lorenz gauge requires the 
following condition on the MP, 



-L k 



0. 



(2.18) 



Using this condition in Eq. (2.171 then provides a wave equation that must be satisfied by the gauge generator 



(2.19) 



A gauge generator that satisfies this equation is unique only up to some 'E!^ that satisfies the homogeneous version 
of (2.19). Specifying the initial data and boundary values (if any) removes the residual gauge freedom and fully 
determines the gauge. 

We consider next the spherical harmonic decomposition of the gauge vector. Momentarily considering again both 
even- and odd-parity, can be broken down into 



im 1 



(2.20) 



There are three even-parity amplitudes and one odd-parity amplitude. We will concern ourselves with determining 
^t, £, r, and in a subsequent paper. In this paper we seek to obtain ^q. Substituting the decomposition of into 
Eq. (2.19), we find after a bit of calculation that satisfies the differential equation 



Wi^o(t,r) = 2/*Rw + /P- 



(2.21) 



Once the gauge generator amplitude is known, we decompose Eq. (2.17) in harmonic amplitudes and see that the 
odd-parity MP amplitudes are transformed by 



(2.22) 







dt ' 






n ~ So; 

or r 


h}^{t,r) 


= -2^0. 





C. Local nature of the metric perturbation and gauge generator at 



The RHS of Eq. (2.21 ) is singu lar at the location of the particle. In this sense Eq. (2.21 ) is very similar to Eq. (|2.10) 
In Ref. [H] we examined Eq. (2.10) to determine the local behavior of V&q. Assuming — 

*o Kit) 



r], we calculated jumps in the field, l^'dp, and in its radial derivative, |9r^olp- We use a subscript p 
to indicate that a function of r is evaluated at the location of the particle, r — rp{t), becoming a function of time. 

Following the same logic here, we postulate a form for the gauge amplitude of £,o = ^[^ ~ ''pCO] ^ ~ A ■ 
Then, similar analysis to that found in Ref. indicates that is C°, i.e. I^olp = 0. Further, we find the jump in 
the first radial derivative is 



fp 



J V 



(2.23) 



Here p{t) comes from the source amplitude P. All three source amplitudes are delta distributions with time dependent 
amplitudes: P = p{t) 5[r - rp(i)], P* = p\t) 5[r - rp(i)], and P'' = p'^(t) 5[r - rp(t)]. 
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Having computed the expected jumps in and its radial derivative, we can use Eq. (2.22 1 to find the jumps in the 
Lorenz gauge MP amphtudes. As with the fields ^'o and ^o, we expect each MP amplitude to consist of left and right 
side differentiable functions that are joined at the location of the particle by Heaviside functions. We calculated the 
jumps m hf" and hf"^ in Ref. The discontinuities in the RW gauge MP amplitudes are exactly canceled out 

by terms arising from the derivatives of and all three Lorenz gauge amplitudes are as expected. The jumps in 
their first derivatives are also analytically computable (either by examining the jumps in the higher-order derivatives 
of or more simply by directly analyzing the Lorenz gauge field equations). We find 

(t) = -jS-r,p\ ldrhr% it) = -J^,V\ \drh^\ (0 = -.j^.V- (2-24) 

J p p J p p J p p 

We use these expressions later (see Fig. [9]) as a powerful check that we have correctly solved the gauge transformation 
equations to high accuracy. 



III. TWO EHS-LIKE METHODS FOR EQUATIONS WITH NON-COMPACT SOURCES 



In Ref. [22] we solved Eq. (2.101 for a variety of eccentric orbits. We used a ED approach to find the Fourier 
harmonic modes of \l/o and transformed back to the TD using the EHS method. The EHS method was first applied 
to wave equations with delta function sources. It allows TD reconstruction of the spherical harmonic amplitudes with 
exponential convergence, circumventing the Gibbs phenomenon that otherwise arises from solving equations with 
discontinuous or singular sources. Our application of the method also demonstrated it could be applied to sources 
with a derivative of a delta function. 

With only a change in spin parameter. 



Eq. (2.21) has a similar differential operator as Eq. (2.101. Where the two 



equations differ markedly is in their source terms. While the source in Eq. (2.10) is point-singular and compact, the 
source in Eq. (2.21 1 is both distributional and non-compact. We can use the linearity of the equation to split off the 
singular part and split the generator into two parts, = Co''* + Co that satisfy separate equations. 



WiCr^(i,r-)=/pp(i)5[r-rp(t)] 



(3.1) 
(3.2) 



While the former equation can be solved using the EHS method, the latter's extended source term is more problematic. 
The extended source is both non-compact and has a time-dependent discontinuity that moves periodically between 
rinin and r,„ax as the particle orbits. In this section we present two equivalent methods for solving Eq. (3.2) using FD 
methods, both of which provide exponential convergence upon returning to the TD. 

As discussed earlier, an eccentric orbit on Schwarzschild provides two fundamental frequencies. When we Fourier 



transform Eq. (3.2 1, we have a two-fold countably infinite frequency spectrum, 



m, n e Z. 



(3.3) 



The Fourier transform and standard TD reconstruction of Co^^{t,r) is then 



T 



(3.4) 



Note that in addition to the already suppressed indices i and m, FD quantities have a third implied index, n. 
Equivalent expressions are used for the Fourier transforms and series representations of the other fields that we 
consider below. Note that while we use the standard tilde notation with Co''* to indicate a FD quantity, for other 
quantities we try to maintain consistency with previous literature by changing the base symbol. For the TD function 
^RWi we write i?Rw in the FD. Similarly, for its TD source term Srw, we write in the FD. 



A. First approach: partial annihilator method 



Our first method for solving Eq. (3.2) is a generalization of the standard method of annihilators used for solving 
inhomogeneous differential equations. It hinges on finding an "annihilator," a differential operator which gives a 
vanishing result after acting on the source. Then, one can act with the annihilator on both sides of the differential 



equation. What results is a homogeneous differential equation of higher order. Our strategy for solving Eq. (3.2) 
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is essentially the same, except that because the initial source is discontinuous the operator that we find does not 
completely annihilate the RHS but instead converts it to a distribution. Hence, we refer to the operator as a partial 
annihilator. 

The RHS of Eq. (3.2) is well suited to the method of partial annihilators because the Regge- Wheeler variable 
satisfies its own wave equation with a point-singular source, (2.151. Therefore, upon dividing Eq. (3.2) by /, we can 
take W2 as the partial annihilator and act on both sides of the equation 



(3.5) 



/W2 ( J Wiat.r) ) = 2fSjm{t,rpit)) = 2/ (GRw(i) S [r - rp{t)] + Fj^it) S' [r - rp{ty 



For simplicity here and in the remainder of this section we drop the ™' tags. We have multiplied back through by / 
to ensure that the leading-order derivatives have unit coefficients. This differential equation is now fourth-order, but 
its source is point-singular. This allow us to solve it using the EHS method, generalized to fourth-order equations. 
The specific form of the source in Eq. (3.5) is given by Martel [32], though we assume that both Grw and Fsw have 
been evaluated at r = rr, 



rp(i). 

We Fourier transform Eq. (3.5) to obtain the FD equation 

1 



fC2 



f 



£ie(r) = 2fZnw{r). 



(3.6) 



The Fourier transform averages the point source motion in time and produces Z^^(r) which has support only within 
the source libration region < r < r^ax- 

There are four linearly independent homogeneous solutions to Eq. (3.6). Two of these are the solutions to the 
cond-order equation £1^ = and we denote ther 
infinity while the other is downgoing at the horizon 



second-order equation £1^ = and we denote them by £,^2- One behaves asymptotically as an outgoing wave at 



i 



h2 



(r^2M), I 



+ 

h2 



(r — )■ Qo). 



(3.7) 



The other two solutions only satisfy the full fourth-order equation £2(7 ^ £iO — 0- such we give them the label 
hA and asymptotic analysis shows that 



(r — ^ 00). 



(3.8) 



These four solutions form a basis spanning the space of homogeneous solutions of Eq. ( |3.6| . The particular solution 
will be a linear combination of these with variable coefficients 



IpM = ^-^(r) 4-2(r) + c+ (r) 1+ (r) + c-Jr) ~^JJt) -(- c+ (r) e+ (r) 



(3.9) 



The four normalization functions c^2/m('') ^'^^ fixed by the method of variation of parameters, which entails solving 
the equations 



dc 



h2/hA 



h2/hA 



(r) 



W{r) 



(3.10) 



Here W{r) is the Wronskian and W^^/hi^^) "modified Wronskian" (Cramer's rule), which is the Wronskian 

with the column corresponding to the C^2//i4(^) homogeneous solution replaced by the column vector (0, 0, 0, 1). Note 
that because the differential operator in Eq. (3.5) is written in terms of r^, the derivatives within the Wronskian must 
also be taken with respect to and the LHS of Eq. (3.10) is a derivative taken with respect to r*. For the two "-f" 

,he fa 



equations, the integral form of Eq. (3.101 is (we change the variable of integration to r and see the factor of / cancel) 

(•Tr 



<2/m(0 - 2 



r Jo 



(GRw(t) S [r' - rp{t)] + FRw(i) S' [r' - rp{t)] 



'dt 



K2/H,y 

W{r') 



-dr'. 



(3.11) 



Likewise, for the two "— " equations (note the change on the limits of integration). 



T, 



r Jo 



(Grw(0 6 [/ - rp{t)] + pRwit) 6' [r' - rp{t)] 



e'^^^dt 



h2/M 



(/) 



W{r') 



dr'. 



(3.12) 
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The EHS method requires knowing only the terminal values of the four functions c^2/?i4('')' ^h2/hi ~ '^?i2/h4(''niax) 
(''mill)- Switching the order of integration and integrating by parts, we find 



T 



Tr 



Guw{t) 



r Jo 



W, 



drW(rp) + 



W{t,) 



At this point we define the EHS in the FD to be 

CM ^ CyJ-2ir) + CyJUr). ^ C+ f+ (r) + C+ 

and the EHS in the TD are defined by the Fourier sums (recall the suppressed £, m, n indices) 



(3.13) 

(3.14) 
(3.15) 



The extension of these solutions to r = rp{t) then gives the desired solution to Eq. (3.2) 

S,T\t. r) = C(t, r)e[r- r^{t)] + C{t, r) 9 [r,{t) - r] . 



(3.16) 



B. Second approach: method of extended particular solutions 



to 



Now we look for a solution to Eq. ( 3.2 ) that does not require a partial annihilator. In the FD the equation transforms 



RW- 



(3.17) 



Again, for notational simplicity we drop the '^^^ tags for the remainder of this section. In the end we want solutions to 
Eq. (3.17) that allow us to form an exponentially converging solution to Eq. (3.2 1 when we transfer to the TD. This 



will require a new technique which we call extended particular solutions, and is closely analogous to the EHS method. 
First, though, we consider how to get the correct causal solution to Eq. (3.171 from a "standard" approach. 

In the subsequent sections we make a distinc tion b etween quantities with °° and ^ tags which designate functions 
computed from a "standard" RHS source (Eq. (3.231 below) and those quantities with + and ~ tags which designate 
functions computed from an "extended" RHS source (Eq. ( 3.32[ ) below). Because the homogeneous solutions do not 
depend on the source, we always tag them with + or ^. We distinguish between particular and homogeneous solutions 
by using the respective subscripts p and /j. 



1. Finding standard FD solutions with causal boundary conditions 

By examining the source, 2/i?ftvvj ^^nd the differential operator, Ci, we can obtain asymptotic and Taylor expansions 
of the particular solution near infinity and the horizon, respectively. The expansions are useful numerically but 
for our purposes here we need only consider the leading asymptotic dependence. (See App.|A]for discussion of the 
asymptotic expansion (r — )■ 00) of and how it couples to the expansion of i?Rw-) 

Consider first the spatial infinity side. Let the RW function have an asymptotic amplitude so i?Rw — ^Ym^^^^* 

as — >■ cxD. Using an asymptotic approximation to 



as — > CXI. We then make the ansatz that = C^re^' 
Eq. (|3.17|) we find 



dr'i 



(C^re"^'-') = 2C] 



RW"= 



lU) 



(3.18) 



Therefore, the asymptotic form of is 



RW' 



r — >■ 00. 



Next we consider the horizon side. The RW function is asymptotically i?Rw 
we expect the particular solution to behave as = Cp Je"*"''* as - 
leading parts of the differential operator we find 



(3.19) 

-00. In this case 



drl 



RW^ 



00. Again, acting with the near-horizon 
/ 1 



M 



RW 



(3.20) 
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Therefore, the near-horizon form of is 



lUJ 

M 



2M. 



(3.21) 



We can use Eqs. (3.19) and (3.21) to set boundary conditions (B.C.'s) for two separate integrations of the inho- 



mogeneous differential equation, (3.17) (yielding two different particular solutions that differ by some homogeneous 



solution). This integration requires the source Rbw: which is itself the solution to the differential equation 



We find it by variation of parameters, which yields 



Kwir) = c+^(r)i?+(r) + c^^{r)R-{r), 



where R^{r) are unit-normalized homogeneous solutions to Eq. (3.22). Note that c^^{r > rmax) 
Crw(^ — ''mill) = Cp^y. Furt hermore, the solution R^^{r) iajwt the same as the EHS to Eq. (3.22). 



Having solved Eq. (3.22) for the source term in Eq. (3.23), and determined the B.C.'s, we are ready to solve 



(3.22) 



(3.23) 



= and 



Eq. (3.17). The idea is to find the two different particular solutions, neither of which has the proper causal behavior. 



and then correct for the acausality by adding appropriate homogeneous solutions. The result of this process is a 
solution to Eq. (3.17) with causal behavior on both sides. The details follow in a series of steps. 



1. Solve for the particular solution from the spatial infinity side (see Fig. [2]). 

We set a B.C. to Eq. (3.17) using Eq. ( |3.19[ ) and integrate to large negative using Eq. (3.23) as the source. 



Although the starting B.C. specified no homogeneous contribution, homogeneous solutions on the horizon side 
will be excited. See the left side of Fig. [2] The particular solution integrated from this side has the asymptotic 
behavior 



-> +00, 

— )■ — oo. 



(3.24) 



The term with coefficient Cp is the part directly dependent on the source and it is sub-dominant in comparison 
to the homogeneous solutions. The coefficients are to-be-determined. Importantly, K+e'"''* is an acausal 
term (upgoing from the past horizon). 




- Real 

- Imaginary 



Outgoing B.C. 



Direction of integration 
< 1 



J I L 



-500 -400 -300 -200 -100 



100 200 300 400 500 600 700 800 900 1000 

r^/M 



FIG. 2. Integration from large r* of the particular solution, Dotted lines indicate the source libration region. 

2. Solve for the particular solution from the horizon side (see Fig. [s]). 



We set a B.C. to Eq. (3.17) using Eq. (3.21) and integrate to large positive using Eq. (3.23) as the source. 



Although the starting B.C. specified no homogeneous contribution, homogeneous solutions on the spatial infinity 
side will be excited. The effect can be seen on the right side of Fig.|3] The particular solution dominates, but the 
constant offset between real and imaginary parts indicates the presence of asymptotically-constant-amplitude 
homogeneous terms. Analysis shows that the particular solution integrated from the horizon will behave as 



-oo, 
-oo. 



(3.25) 
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The coefficients are to-be-determined and again we find an acausal term (A e which in this case is 

ingoing from past null infinity. 




200 300 



900 1000 



FIG. 3. Integration from large negative r» of the particular solution, . Dotted lines indicate the source libration region. 
3. Solve for the homogeneous solution from the spatial infinity side. 



Next, we set an outgoing B.C. to the homogeneous version of Eq. (3.171 at large positive r*. We integrate to 



solve the scattering problem for reflection and transmission amplitudes and 33]. In the terminology of 
Gal'tsov i34J this is an "up" mode. 



+ „iu>r. 



T+e 



R+e 



-> -l-oo, 
— oo. 



(3.26) 



Scaled appropriately, this solution can be added to Eq. (3.24 1 to remove its acausality. 



4. Solve for the homogeneous solution from the horizon side. 



the scattering problem for reflection and transmission amplitudes R and T . This is an "in" mode 



We set a downgoing B.C. to the homogeneous version of Eq. ( |3.17[ ) at large negative and integrate to solve 

(3.27) 



T-e- 
R-eJ' 



-oo, 

-CO. 



Scaled appropriately, this solution can be added to Eq. ( 3.25[ ) to remove its acausality. 
5. Resolve the acausality in the particular solutions (see Fig.[4]). 



The acausal piece in Eq. (3.24 1 is K+e 



Using Eq. (3.261, we can remove this by subtracting k^^, 



+ K e 



K+R+t 



-oo. 



(3.28) 



The acausal piece in Eq. (3.251 is A e "^''*. Using Eq. (3.27), we can remove this by subtracting A ^ 



A-T-e^'"'^*, 



— — oo, 
— ^ -foo. 



(3.29) 



Eqs. (3.28) and ( 3.29[ ) are both solutions to Eq. (3.17) and both satisfy the causal nature of the problem. 
Therefore they must be equal. In order to form them, we must know k+ and A~. We find them by setting 
Eqs. (3.28) and (3.29) and their first derivatives equal at any point. 



->^~Ch (3-30) 

a.if - = drj^ - K+drJ+. (3.31) 

We solve these equations for k+ and A^ and form — A^^^ and — k'^^,^ , which are equivalent. 



The function ^^'^ — — ^ £.h ~ ^ ^^it represents the standard solution to Eq. (3.171. If the TD source 
were differentiable, we would be able to find the corresponding TD solution via an exponentially converging Fourier 
synthesis. However, the source in this case is non-differentiable and we need an EHS-like trick to complete the method. 
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FIG. 4. Causally correct solution to Eq. (3.171, ^ . Dotted lines indicate the source libration region. 



2. Restoring exponential convergence with extended particular solutions 



The EHS of the Regge- Wheeler equation Eq. (3.22) are found by taking the constants and scaling the unit- 



normalized homogeneous solutions 



i?±(r). 



(3.32) 



These solutions are defined for all r > 2M. In like fashion we seek to find FD EPS of Eq. (3.171 and denote these 



by We first find by separately integrating Eq. (3.171 with the modified source terms -R^^. The solutions are 
each made to match the exterior behavior of by adding the correctly scaled homogeneous solutions found in Step 
[5] above. We then define 



i+ EE 4+ - t ^ - A-c 

See Fig.[5j which contrasts Fig. |4]in the source region. 

These FD EPS can be transferred to the TD via Fourier series 



(3.33) 



(3.34) 



The solution to Eq. (3.2) is then the weak solution, 

r) = e+(t, r)0[r~ rp{t)] + Cit, r) 6 [r,it) - r] 



(3.35) 



The support for this claim has three legs. Firstly, the same arguments about EHS, based on analytic continuation, 
made by Barack, Ori, and Sago in Ref. |19| appear to apply in extension to Eq. (3.2) as well. Secondly, we demonstrate 



existence numerically by integrating the equation, with causal boundary conditions, and checking that the jump 
conditions (internal boundary conditions) at the particle are satisfied. One then appeals to the linearity of the 
equation to establish uniqueness. Finally, we have an independent numerical solution found through the method 
of partial annihilators and given in Eq. ( |3.16 1. We have confirmed that the two methods give entirely consistent 
solutions. These results are covered in detail in the next section. 



IV. RESULTS 



The methods of the previous section allow us to transform odd-parity solutions of the first-order Einstein equations 
from RW to Lorenz gauge. As an example, we consider an orbit with eccentricity e — 0.764124 and semi-latus rectum 
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_2 5 I I I I I ^ I I I ' I ' I ' I ' I ' I > I > I 

-500 -400 -300 -200 -100 100 200 300 400 500 600 700 800 900 1000 



FIG. 5. Causally correct extended particular solution, Note the difference within the libration region (shown in dotted 
lines), between and (^^^'^ (Fig. Ul. The difference in TD convergence between these two solutions is shown later in Fig. Isl 




r^/M 

FIG. 6. The £ = 2,m — 1 mode of the odd-parity RW-to-Lorenz gauge generator amplitude £,o- The Lorenz gauge MP 
amplitude /i2 differs from only by a factor of —2. Dotted lines indicate the region of libration. Orbital parameters are given 
in the text. While the field /i2 grows asymptotically, upon transforming to an orthonormal frame it would contribute a term 
that falls off as 1/r. 



p = 8.75455. This orbit was used in Fig. [T] where we showed the RW amplitudes ht and for £ = 2 and m = 1 
(/i2 = 0). The MPs can be evaluated at any time but we chose to display results at t = 93.58 (where t = is at the 
periapsis). The RW modes are discontinuous at r = rp{t) and lack asymptotic flatness. The gauge generator to go 
from RW to Lorenz gauge is co mpute d for this same orbit and at the same time in the TD. It is used to obtain the 
MPs in Lorenz gauge using Eq. ( 2.22[ ). Fig. [6] shows the £ = 2, m — 1 amplitude of the gauge generator itself, which 
differs from h2 in Lorenz gauge only by a factor of —2. Fig. [t] shows the Lorenz gauge metric amplitudes h]' and h]; 
for the same mode. The MPs are now C*' at r = rp{t) and are asymptotically flat. 

Of key importance to our method is the exponential convergence of the TD solutions. We can first consider self- 
convergence of the modes for all r. As an example we choose an orbit with e = 0.188917 and p = 7.50478 at time 
t = 96.44. In Fig. ID we show the self-convergence of £,o{r) for a set of partial Fourier sums over n from —N < n < +N 
for various N. The right panel of this flgure shows exponential self-convergence of the EPS method as a function of 
r, including at the particle. This result is in contrast with the left panel which shows that the standard method is 
only algebraically convergent in the source libration region. Note that the convergence is initially exponential before 
becoming algebraic around an error level of 10~*. This transition is due to the equation for having singular and 
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FIG. 7. The £ — 2,m = 1 mode of the Lorenz gauge MP amphtudes and h^. Dotted hnes indicate the region of libration. 
Orbital parameters are given in the text. Note (comparing to Fig. [TJ the discontinuity at the location of the particle has 
vanished and the wave no longer grows asymptotically. We plot fh^ so we can see the wave behavior near the horizon. 



extended parts (see Eqs. (3.1) and (3.2)). We find the singular part using EHS, which converges exponentially. This 
part of the solution dominates the self-convergence in the left panel at first. Eventually, the lack of differentiability 
of the extended source and the use (for comparison) of the standard Fourier series for that part of the field manifests 
itself. The appearance of Gibbs behavior stalls the convergence in the libration region. 

Beyond self-convergence, we can check absolute convergence to the analytically-known jump conditions. This test 
can be applied to both the MP amplitudes and their first radial derivatives. As shown in Fig. [9j we find exponential 
convergence to the analytically computed values in Eq. (2.24), in this case using the partial annihilator method. In 
the left panel we consider a less eccentric orbit (e = 0.188917 and p = 7.50478) and in the right panel a more eccentric 
one (e = 0.764124 and p — 8.75455). Each partial Fourier sum ranges over all harmonics from —N < n < +N for 
different values of N. We also compute the partial sums at several points in time during an orbit and compare the 
discrepancy between the expected jump and the numerically computed jump. Since the Lorenz gauge amplitudes are 
all C°, we plot absolute convergence for the jumps in the amplitudes themselves (which are expected to converge to 
0) and relative convergence for the jumps in the r derivatives of the amplitudes. As expected the low eccentricity 
orbit converges with fewer terms and yet there is no difficulty reaching convergence in the high eccentricity case as 
well. In each case the convergence appears to bottom out around 10~^^. 

Although we have only displayed results in this paper for the £ = 2, m = 1 mode, we have run the code on many 
different modes and for different orbits. We have no difficulty in computing the gauge transformation from RW to 
Lorenz for odd-parity modes with high accuracy. It now remains for us to apply these methods to the even-parity 
part of the gauge transformation, a somewhat more involved procedure. We will turn to that issue in a subsequent 
paper. 



V. CONCLUSION 



This paper is the first of two on the transformation of metric perturbations from Regge- Wheeler gauge to Lorenz 
gauge. This first paper was confined to treating the odd-parity part of the MPs and devoted much of the discussion 
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FIG. 8. Self convergence of the £ = 2, m = 1 mode of ^o- We show results for the standard method in the left panel and 
our (equivalent) partial annihilators and EPS on the right. The orbital parameters are e — 0.188917 and p = 7.50478, and the 
fields are calculated at time t = 96.44. 




FIG. 9. Convergence of the jump conditions in the values of the MP amplitudes and their first radial derivatives. In the left 
panel we show results for an orbit with e = 0.188917 and p = 7.50478. On the right we have e — 0.764124 and p — 8.75455. 



to the development of two new analytic/numerical methods for using frequency domain methods to find accurate 
sohitions in the time domain. The follow-on paper will be primarily devoted to discussing the analytic problem of 
finding the even-parity part of the gauge transformation, and will draw upon the numerical methods which we have 
detailed here. 
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Appendix A: Asymptotic expansions and boundary conditions 

In the RWZ formahsm it is useful to compute asymptotic expansions of the master functions about r — co to provide 
boundary conditions for starting numerical integrations at finite radius. In this paper, the inhomogeneous equation 
for the gauge generator, Eq. (3.17), has a source term that is non-compact. This fact leads to an inhomogeneous 



recurrence relation for the asymptotic expansion of that requires as input the asymptotic expansion of the source 
term. 

We start by writing 



(Al) 



where Jo{r) is the Jost function [33], which goes to 1 at infinity. We use Eq. (2.14) to express the RHS of Eq. (3.17) 



in terms of the CPM function. Then we Fourier transform that function and plug in Eq. (Al) to obtain 

£{£- 



2 1 



M 
r 



dr 



Jo 



2iuj 



2M 



1) 



Jo — —i^JR- 



Here Jr — J^^„„ from App. D of 22J. Now, we assume the following forms of J^, and J^, 



(A2) 



3=0 ^ ' 3=0 ^ ' 

Plugging these in and assuming the equation is satisfied order-by-order gives the coupled recurrence formula. 



2z(j-l)a°= (j-2){j~ I) -£{£+!) 



a°_i + 2Mu 



(A3) 



(A4) 



ttj , given in Eq. D5 of | 



Assuming a° = for j < 0, this recurrence allows for the calculation 
1, which represents the homogeneous solution to Eq. (A2). We can 



The coefficients 

of all a°. Note that this recurrence fails at j 
choose that coefficient to be anything. ' ' 

The particular solution here is identical to the homogeneous solution to the fourth-order equation, given asymptot- 
ically in Eq. (3.8). We can use this asymptotic expansion for both situations. 



On the horizon side, where the potential falls away exponentially, it is enough to use the expression in Eq. (3.20) 
and a sufficiently large and negative starting location for integration. A Taylor expansion could be used if the 
starting location were farther from the horizon. The boundary conditions to the second-order homogeneous solutions 
are exactly analogous to those given in the odd-parity recurrence of App. D in [2?. The only difference is a change 
of the spin parameter in the potential from 2 to 1. 
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